#################  load R packages --- needs to stay here 
rm(list=ls())   #### cleans memory - needs to stay here
#args <- commandArgs(trailingOnly = TRUE)
############  Gridcell Simulation Switch  ##########################################################
GridsimulationSwitch=c('OFF','ON')[2]  
########1=single point simulation, 2= Grid cell simulation
start_time <- Sys.time()

###########load packages######################################################################
list.of.packages <- c("ggplot2", "plyr","parallel")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(ggplot2)
library(plyr)
library(parallel)
library(sp)
# Check if the 'SPEI' and 'PDSI' packages are installed
if (!requireNamespace("SPEI", quietly = TRUE) | !requireNamespace("scPDSI", quietly = TRUE)) {
  warning("The 'SPEI' or 'PDSI' package is not installed. Please install it to proceed.")
} else {
  library(SPEI)
  library(scPDSI)
}

if(dev.cur()>1){replicate(dev.cur()-1,dev.off())} 
workDir <- c("/discover/nobackup/myang5/VACS/SIMPLE_DISCOVER/")
#setwd(script.dir <- dirname(sys.frame(1)$ofile)) #### Relative working directory
setwd(workDir)  ########Absolute working directory

climate_period_list <- c("baseline", "historical", "future") # Add your list of climate periods
#climate_period_list <- c("historical")
MapExtention_list <- c("Algeria","Angola","Benin","Botswana","Burkina Faso",                                    # 5
                       "Burundi","Cameroon","Central African Republic", "Chad","Comoros",                       #10
                       "Democratic Republic of the Congo","Djibouti","Egypt", "Equatorial Guinea","Eritrea",    #15
                       "Ethiopia","Gabon","Gambia", "Ghana","Guinea",                                           #20
                       "Guinea Bissau","Ivory Coast","Kenya", "Lesotho","Liberia",                              #25
                       "Libya","Madagascar","Malawi", "Mali","Mauritania",                                      #30
                       "Morocco","Mozambique","Namibia",  "Niger","Nigeria",                                    #35
                       "Republic of the Congo","Rwanda","S. Sudan", "Sao Tome and Principe","Senegal",          #40
                       "Sierra Leone","Somalia","Somaliland", "South Africa","Sudan",                           #45
                       "Swaziland","Togo","Tunisia",  "Uganda","United Republic of Tanzania",                   #50
                       "Western Sahara","Zambia","Zimbabwe")
#MapExtention_list <- c("South Africa")

#species_list <- c("fonio", "tef", "sorghum", "soybean", "groundnut", "sesame", "taro", "tomato") # Add your list of species
#species_list <- c("fingermillet","soybean", "cowpea","grasspea","pigeonpea","lablab","bambaragroundnut","cassava","yams","okra","africaneggplant","josephscoat")
#species_list <- c("fonio", "tef", "sorghum", "fingermillet", "okra", "josephscoat")
#species_list <- c("groundnut", "maize", "sesame", "tomato", "taro", "cocoyam", "sweetpotato", "cassava")
#species_list <- c("sweetpotato", "soybean", "pigeonpea", "maize", "lablab", "groundnut", "grasspea", "cassava", "cowpea")
#species_list <- c("maize","fonio","tef","fingermillet","sorghum","pearlmillet","plantain","pumpkin","soybean","cowpea","grasspea","mungbean","pigeonpea","lablab","bambaragroundnut","groundnut","sesame","cassava","cocoyam","taro","sweetpotato","yams","tomato","okra","africaneggplant","josephscoat")
species_list <- c("fonio","tef","fingermillet","sorghum","pearlmillet","plantain","pumpkin","soybean","grasspea","mungbean","pigeonpea","lablab","bambaragroundnut","groundnut","sesame","cassava","cocoyam","taro","sweetpotato","yams","tomato","okra","africaneggplant","josephscoat")
#species_list <- c("maize","fonio","sorghum","soybean","cowpea","bambaragroundnut","cassava","cocoyam","sweetpotato","yams")
cultivar <- 1
WeatherType=c("WTH","CSV","Rdata","AgMIP")[4]

for (species in species_list) {
for (climate_period in climate_period_list) {
  
  if (climate_period == "baseline") {
    climate_scenario_list <- c("ERA")
    climate_model_list <- c("5")
  } else if (climate_period == "historical") {
    climate_scenario_list <- c("SSP370")
    climate_model_list <- c("GFDL-ESM4", "IPSL-CM6A-LR", "MPI-ESM1-2-HR", "MRI-ESM2-0")
  } else {
    climate_scenario_list <- c("SSP126","SSP370")
    climate_model_list <- c("GFDL-ESM4", "IPSL-CM6A-LR", "MPI-ESM1-2-HR", "MRI-ESM2-0")
  }

  for (climate_scenario in climate_scenario_list) {
    for (climate_model in climate_model_list) {
      for (MapExtention in MapExtention_list) {
          setwd(workDir)
          # Check if the output file exist or not. If not, create a dummy file.
          outFileName <- paste0("./Output/Gridcell_summary_",species,"_",climate_period,"_",climate_scenario,"_",climate_model,"_",MapExtention,".csv")
          if(file.exists(outFileName)) {next} else {
            dummy_df <- data.frame(matrix(vector(mode = "numeric", length = 0L), ncol = 14))
            colnames(dummy_df) <- c("Exp.","Species.","Trt.","row","col","long","lat","SowingDate","Duration","Biomass","Yield","MaturityDay","SowingYear","MaturityYear")
            write.csv(dummy_df,outFileName,row.names = F)
          }
          
          print(outFileName)
          source("Mainfunction.R")
          
          if(climate_period == "baseline"){
            WeatherDir=paste0("./Gridcell Weather/",climate_period,"/")
          } else if (climate_period == "historical") {
            #climate_scenario<-"SSP370"  #only used for renaming historical output file because 2015-2019 was filled with SSP370
            WeatherDir=paste0("./Gridcell Weather/",climate_period,"/",climate_model,"/")
          } else {
            WeatherDir=paste0("./Gridcell Weather/",climate_period,"/",climate_scenario,"/",climate_model,"/")
          }
          
          ###set years to simulate###
          if(climate_period == "baseline"){
            SimulatingYear=c(2000:2019)
          } else if (climate_period == "historical") {
            SimulatingYear=c(1990:2019)
          } else {
            SimulatingYear=c(2035:2064)
          }
          #to test smaller subset of years, uncomment out next line
          #SimulatingYear=c(2010)
          
          #########AutoGenGridInFile##########################
          if (GridsimulationSwitch == 'ON') {
            source("AutoGenGridInFile.R")
            GridInData <- AutoGenGridInFile(workDir, WeatherDir, MapExtention, species, cultivar)
            if (prod(dim(GridInData)) == 0) next
          }
          
          ########## Output option################
          DailyOutputforgridcell=c('OFF','ON')[1]
          DailyOutputOutput=c("Crop","Exp","Label","Trt","Day","DATE","Tmax","Tmin","Radiation",
                              "TT","fSolar","Biomass","dBiomass","HI","Yield","F_Temp","F_Heat",
                              "F_Water","ARID","I50B","I50A","ETO","MaturityDay")
          
          #if(climate_period == "baseline"){
          #  MapoutputYear=2010
          #} else if (climate_period == "historical") {
          #  MapoutputYear=2010
          #} else {
          #  MapoutputYear=2050
          #}
          #to plot different year, uncomment out next line
          #MapoutputYear=2010
          
          
          ############  Model starts here  ##################################################################
          ############  Read inputs - treatments, soil weather, CO2
          if(GridsimulationSwitch=='OFF'){
            management<-read.table("./Input/Simulation Management.csv",header=TRUE,sep=",",stringsAsFactors=FALSE)
            treatment<-read.table("./Input/Treatment.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);treatment$Species.=tolower(treatment$Species.)
            cultivar<-read.table("./Input/Cultivar.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);cultivar$Species.=tolower(cultivar$Species.)
            irri<-read.table("./Input/Irrigation.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);irri$Species.=tolower(irri$Species.)
            soil<-read.table("./Input/Soil.csv",header=TRUE,sep=",",stringsAsFactors=FALSE)  
            para_spec<-read.table("./Input/Species parameter.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);para_spec$Species.=tolower(para_spec$Species.)
            
            ####match experiment
            management<-management[management$ON_Off==1,]
            treatment<-merge(treatment,management,by=c("Species.","Exp.","Trt."), suffixes=c("",".y"))
            treatment<-merge(treatment,soil,by="SoilName.")
            treatment<-merge(treatment,para_spec,by="Species.")
            treatment<-merge(treatment,cultivar,by=c("Species.","Cultivar."))
          }else{
            setwd(workDir)
            #treatment<-read.table("./Input/Grid_input_Ghana_baseline.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);treatment$Species.=tolower(treatment$Species.)
            #treatment<-read.table("./Input/Grid_input_Ghana_future_ssp126_gfdl-esm4.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);treatment$Species.=tolower(treatment$Species.)
            #treatment<-read.table(GridInFileName,header=TRUE,sep=",",stringsAsFactors=FALSE);treatment$Species.=tolower(treatment$Species.)
            treatment<-GridInData;treatment$Species.=tolower(treatment$Species.)
            para_spec<-read.table("./Input/Species parameter.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);para_spec$Species.=tolower(para_spec$Species.)
            irri<-read.table("./Input/Irrigation.csv",header=TRUE,sep=",",stringsAsFactors=FALSE);irri$Species.=tolower(irri$Species.)
            treatment<-merge(treatment,para_spec,by="Species.")
            treatmentsingle<-treatment
            x=1:length(SimulatingYear)
            no_cores <- detectCores()
            cl <- makeCluster(mc <- getOption("cl.cores", no_cores))
            clusterExport(cl, c("treatment","SimulatingYear"))
            results <- parLapply(cl,x,Treatmentplusyear) 
            stopCluster(cl)
            treatment <- do.call('rbind',results)
          }
          
          
          
          
          
          
          RunModel=function(i){
            source("Mainfunction.R")
            paras=ParaInput(i)
            res<-tryCatch({SIMPLE(para=paras[c(1:3)],weather=paras$weather,ARID=paras$ARID)},error=function(e){cat("ERROR :",conditionMessage(e),"\n")})
            return(res)
          }
          
          
          
          
          
          
          #####Run all the experiment one by one
          # x=1:nrow(treatment)
          # results=list()
          # if(GridsimulationSwitch=='OFF'){observations=list()}
          # for (i in 1:length(x))
          # {
          #   results[[i]]=list()
          #   source("Mainfunction.R")
          #   res=RunModel(x[i])
          #   results[[i]]<-res
          # 
          #   if(GridsimulationSwitch=='OFF'){
          #     source("Obsfunction.R")
          #     obs=ObsInput(x[i])
          #     observations[[i]]<-obs}
          # }
          
          ##########
          
          
          
          
          
          ########parallel running
          t1=Sys.time() 
          x=1:nrow(treatment)
          no_cores <- detectCores()
          cl <- makeCluster(mc <- getOption("cl.cores", no_cores))
          clusterExport(cl, c("treatment","irri","GridsimulationSwitch","WeatherType","WeatherDir","DailyOutputOutput"))
          results <- parLapply(cl,x,RunModel) 
          if(GridsimulationSwitch=='OFF')
          {
            source("Obsfunction.R")
            observations<- parLapply(cl,x,ObsInput) 
          }
          
          stopCluster(cl)
          Sys.time()-t1
          ###########
          
          end_time <- Sys.time()
          print(end_time - start_time)
          
          #########Simulation results reorganization
          res.df <- do.call('rbind',results) 
          Res_daily=ldply(res.df[,1])
          Res_summary=ldply(res.df[,2])
          
          
          if(GridsimulationSwitch=='OFF'){
            obs.df=do.call('rbind',observations)
            Obs_Biomass=ldply(obs.df[,1])
            Obs_FSolar=ldply(obs.df[,2])
            Obs_Yield=ldply(obs.df[,3])
            Obs_Summary=ldply(obs.df[,4])
            Obs_Summary=Obs_Summary[,c('Trt','Exp','Yield')]
            Res_Summary=merge(Res_summary,Obs_Summary,by=c('Trt','Exp'))[,c("Crop","Exp","Label","Yield.x","Yield.y","Duration","Trt")]
            names(Res_Summary)[4:5]=c("Sim_Yield","Obs_Yield")
            ##### write the simulations into files
            Yeargap=paste(unique(format(Res_summary$MaturityDay,"%Y")),collapse = "_")
            Speciesgap=paste(unique(Res_summary$Crop),collapse = "_")
            write.table(Res_daily,paste("./Output/Res_daily_",Speciesgap,"_",Yeargap,".csv",sep=""),col.names=TRUE,row.name=FALSE,sep=",")
            write.table(Res_summary,paste("./Output/Res_summary_",Speciesgap,"_",Yeargap,".csv",sep=""),col.names=TRUE,row.name=FALSE,sep=",")
            
            source("Plot.R")
            gplot(Res_daily,Res_Summary,Obs_Biomass,Obs_FSolar)
          }else{
            treatmentsub=treatmentsingle[,c('Exp.','Species.','Trt.','row','col','long','lat')]
            Res_summary=Res_summary[,c("Exp","SowingDate","Duration","Biomass","Yield","MaturityDay","HeatDays","ARIDDays")]
            Res_Summary=merge(treatmentsub,Res_summary,by.x="Exp.",by.y = "Exp")
            uniqueYears = as.integer(unique(format(Res_Summary$MaturityDay[c(1:length(Res_Summary$MaturityDay))],"%Y")))
            Yeargap=paste(as.character(c(min(uniqueYears),max(uniqueYears))),collapse = "_")
            #Yeargap=paste(unique(format(Res_Summary$MaturityDay[c(1,length(Res_Summary$MaturityDay))],"%Y")),collapse = "_")
            if(DailyOutputforgridcell=='ON'){
              write.table(Res_daily,paste("./Output/Gridcell_daily_",species,"_",climate_period,"_",climate_scenario,"_",climate_model,"_",MapExtention,".csv",sep=""),col.names=TRUE,row.name=FALSE,sep=",")
            }
            Res_Summary$SowingYear=substr(Res_Summary$SowingDate,1,4)
            Res_Summary$MaturityYear=substr(Res_Summary$MaturityDay,1,4)
            setwd(workDir)
            write.csv(Res_Summary,outFileName,row.names = F)
            
          #  Res_SummaryMap=Res_Summary[Res_Summary$SowingYear==MapoutputYear,]
              
          #  source("MapPlot.R")
            
            
          #  Worldcountry=readRDS(file = "./Input/Map/Worldcountry.rds")
          #  Worldstate=readRDS(file = "./Input/Map/ne_10m_admin_1_states_provinces.rds")
          
          #  if(any(MapExtention=="World" | MapExtention=="Africa" )){
          #    WMap <- Worldcountry
          #  }else{WMap<-Worldstate[Worldstate@data$admin%in%MapExtention,]
          #  WMap<- subset(WMap,!fips %in% c("US02", "US15", "US72"))}
            
          #  MapPlot(Res_SummaryMap,index="SowingDate",WMap,MapExtention,Title=paste0(Res_Summary$Species.[1]," sowing date"),Unit=" ")
          #  MapPlot(Res_SummaryMap,index="Duration",WMap,MapExtention,Title=paste0(Res_Summary$Species.[1]," duration"),Unit="Season (days)")
          #  MapPlot(Res_SummaryMap,index="Biomass",WMap,MapExtention,Title=paste0(Res_Summary$Species.[1]," biomass"),Unit="Biomass (kg/ha)")
          #  MapPlot(Res_SummaryMap,index="Yield",WMap,MapExtention,Title=paste0(Res_Summary$Species.[1]," yield"),Unit="Yield (kg/ha)")
          }
        }
      }
    }
  }
}

